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Abstract. A linear algebraic solution is provided tor the problem of retrieving the 
location and time of occurrence of lightning ground strikes from an Advanced 
Lightning Direction Finder (ALDF) network. The ALDF network measures field 
strength, magnetic bearing, and arrival time ot lightning radio emissions and 
solutions for the plane (i.e.. no Earth curvature) are provided that implement all ot 
these measurements. The accuracy of the retrieval method is tested using computer- 
simulated data sets and the relative influence ot bearing and arrival time data on the 
outcome of the final solution is formally demonstrated. The algorithm is sufficiently 
accurate to validate NASA s Optical Transient Detector (Ol D) and Lightning 
Imaging System (LIS). We also introduce a quadratic planar solution that is useful 
when onlv three arrival time measurements are available. The algebra of the 
quadratic root results are examined in detail to clarify what portions ol the analysis 
region lead to fundamental ambiguities in source location. Complex root results are 
shown to be associated w ith the presence of measurement errors when the lightning 
source lies near an outer sensor baseline ot the ALDF network. For arbitrary 
noncollinear network geometries and in the absence of measurement errors, it is 
shown that the two quadratic roots are equivalent (no source location ambiguity) on 
the outer sensor baselines. The accuracy of the quadratic planai method is tested 
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with computer-generated data sets and the results ate generally better than those 
obtained from the three-station linear planar method when bearing errors are about 
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1. Introduction 

Advanced Lightning Direction Finder (ALDF) sensors, developed by Global Atmospherics 
Inc. (GAI), have the ability to detect the Field strength, magnetic bearing, and arrival time of 
lightning radio emissions. In 1992, Lightning Location and Protection, Inc. (a division of GAI) 
completed development of an IMproved PcrformAncc from Combined Technology (IMPACT) 
method for determining the location and time of occurrence of lightning return strokes from these 
data (Cummins et al„ 1993). The IMPACT algorithm is based on minimizing a z 2 function similar 
to that provided in (1) of Hiscox et al. (1984) but generalized to accommodate arrival time data. 
The lightning time of occurrence, /, and the longitude, A, and the latitude, cp, ot the lightning 
source on an ellipsoidal Earth is estimated. 

Since the IMPACT algorithm uses a numerical approach to determine the absolute 
minimum of the nonlinear Z 2 hypersurface, it does not represent an analytic solution to the 
problem, i.e., the source location and time of occurrence are not directly determined in terms of the 
measurements and measuring network geometry. Instead, the algorithm uses the computer to 
search for the optimum values of (A, (p, t) that minimize each new set of measurements implies 
starting an entirely new search tor an answer. In effect, the actual solution is estimated using the 

constraints of the data and the power of the computer. 

Relative multiple minima in the / (A.<p.r) hypersurface can lead to a premature termination 
in the computer search and an erroneous solution. The presence of data errors can generate 
additional relative minima and additional solution errors. Since it is computationally expensive to 
check for erroneous solutions during real-time processing of ALDF data, lightning source 
solutions depend, in general, on where in the solution space (A. cp, t) the computer search begins. 

Although the IMPACT algorithm has had successful practical application (Cummings 
et al., 1998), the specific algorithm software is proprietary and is not widely distributed free of 
charge to the scientific community. In this study, our interest was to develop an economical, yet 
useful (four-station) ALDF ground-truth site at the Melville Island/Darwin, Australia region, and 
elsewhere in the world, that could be used to validate NASA’s space-based lighting detectors (the 
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Lightning Imaging Sensor (LIS) and the Optical Transient Detector (OTD) described in Christian 
ct al. (1992) and Goodman et al. (1995), respectively). We desired a computationally quick 
algorithm that easily ingests all types of ALDF measurements, and that produces ground strike 
locations with accuracies better than 4 km (the nadir resolution of LIS) for sources that are within a 
few hundred kilometers of the ALDF network. This has been the primary motivation behind our 
algorithm development. As such, we did not require an algorithm that accounts for Earth curvature. 

A secondary motivation behind this work was realized during our algorithm development. 
We learned that a serious pedagogical effort to examine the mathematical foundations of lightning 
ground-strike retrievals does not appear in the literature. Perhaps the wide use and practicality of 
powerful numerical computer algorithms have offset the desire to obtain more exact analytic 
solutions. In retrospect, the proper mathematical approach is to first solve the least difficult form of 
the retrieval problem in a thorough way, and then methodically upgrade the formalism to account 
for additional physical complexities such as associated with the addition of new data sets. Earth 
curvature, and wave propagation effects. An effort must also be made to clearly explain 
fundamental ambiguities that arise in the retrieval process when an experimenter has limited 
measurements available due to sensor failure or other causes. Although there have been some 
analytic efforts to account for Earth sphericity for magnetic bearing data inversions (e.g., see 
Orville, 1987) a complete pedagogy that ultimately leads to rigorous retrievals of ground-strikes on 
an ellipsoidal Earth surface (and that uses all ALDF data constraints in a coherent/simultaneous 
and optimum fashion) is notably absent in the literature. 

To the best of the authors’ knowledge, no one has yet explicitly provided an analytic 
solution to the problem of determining lightning source location and time of occurrence using 
collective measurements of field, magnetic bearing, and arrival time measurements when Earth 
curvature and propagation effects are neglected. This is a fundamental starting point. For multiple 
data set inversions, the relative importance of bearing and arrival time data on the outcome of the 
final solution has not been formally demonstrated. Furthermore, no thorough investigation of 
solution ambiguities have been provided when one is limited to just three arrival time 



measurements. This writing introduces theoretical derivations that address each of these problems. 
In so doing, we arrive at the required OTD/LIS ground truth algorithms. 

We specifically determine the source location (.v,y) as a mathematical function of the 
measurements under ot variety of conditions (i.e., ditlcrcnces in the number, location, and type of 
measurements). Since algebraic solutions are obtained, we do not need to invoke a computer 
search algorithm to determine optimum solution parameters. We provide unique physical insight 
into the nature of the retrieval problem because we determine exactly how the measurements arc 
specifically related to the lightning source location (and lime ol occurrence). 

A Linear Planar (LP) method is first introduced that allows one to simultaneously analyze 
field, bearing, and arrival lime measurements. To the best of our knowledge, this has not 
previously been shown in the literature. The method involves one large system of linear equations 
that offers a high degree of flexibility from the point of view of the user’s applicational needs. For 
example, if only a certain number and type of measurements are available in an experiment, the 
linear system of equations degenerates into a smaller set of equations, and a straightforward 
solution process is retained. 

We also introduce a Quadratic Planar (QP) method that can be used when only three arrival 
time measurements are available. Such a situation arises if there are sensor hardware failures and/or 
when field amplitude and bearing measurement data quality is unacceptable. Although this method 
is mathematically nonlinear, full analytic solutions arc derived. Physical insight about the nonlinear 
solution space, not discernible from conventional % 2 analyses, is fully described by examining in 
detail all quadratic root solutions derived from the QP method. For example, we show explicitly 
that a certain mathematical discriminant vanishes for certain lightning locations, and that these 
source locations produce complex roots (negative discriminants) in the presence of measurement 
errors. 

Extensive tests of the LP and QP retrieval methods are provided using computer-simulated 
data sets and these methods are applied in a study of ALDF data that were obtained from the 
Maritime Continent Thunderstorm Experiment (MCTEX) analysis region in Darwin, Australia 
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(Kenncn ct al„ 1994; 1996). Data from this network comprise one of several ground-truth sites for 
the validation of OTD and LIS. 


2. Linear Planar (LP) Method 

We begin by considering 3 sensors situated at locations r, , i — 1,2, ... , /i relative to 
some origin. Each sensor has the capability to measure the arrival time, t , , magnetic bearing, d, , 
and field strength, F, , of the radio emissions from a lightning source with location, r, time of 

occurrence, t, and radiation source strength, s. Hence, fiom the 3 u measurements 
{(/,,$ ,F t ),. <P„,F„) ) we wish to determine the five unknowns (r, /, .v). In so doing, we 

neglect Earth curvature. 

Figure 1 summarizes the geometry of the LP model. Because ALDF sensors might not be 
deployed on a flat topography, the < lh sensor located at r, need not lie in the .tv plane, i.e., c, * 0 in 
general. The relative position vector follows standard physics convention, that is, it points from the 
source at r to the observation point, r,, so that R, = r, - r. Neglecting refractive effects in the 
atmosphere, the excitation time of the i th sensor is 

/.=/ + -/? , (» 

1 c 1 

where c is the speed of light. Solving for . squaring, and rearranging terms leads to 



- c 2 r 2 ) = x r x + y,y + 


- c 2 tJ - -^(r 2 -c 2 t 2 ) . 


( 2 ) 


It is desirable to remove the last term on the right-hand side of (2) since it is nonlinear in the space 
and time variables. To do this, we define the measurement 


a^Ur 2 -c 2 t 2 )-Ur 2 -c 2 t 2 ) 


(3) 
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A comparison of (2) and (3) shows that a, is linearly related to the lightning location, 
r = (x.y,z), and lightning time of occurrence, t, that is, 


a = (.t - .v, )x + (y. - y, )y + (z i -z l )z- c 2 ( r - /, )t ; i = 2 , 3 n 


(4) 


A detailed investigation of this linear form has been provided in Koshak and Solakiewicz (1996). 

Next, we consider the information content of ALDF bearing data. From Figure 1 we see 
that the lightning location (x, y) is given by 


.v = „Y ( + p, cosp, 
v = v, + p, sin <p t , 


(5) 


where p, is the horizontal distance from the ( lh site to the lightning ground strike location. 


It is useful to define the measurement 


/} = X' sinp, - y, cosfy . 


(6) 


Now from (5) we note that (x~x,)/(y - y,) = cos0,/sin</>, , so that using (6) gives 


ft =sin0,x - cosfyy . 


(7) 
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Finally, we consider measurements of the radiated field strength. Assuming a I//?, 
attenuation in the radiation field gives 


F - = T 


( 8 ) 


Once again, we solve for R, , square, rearrange terms, and define the measurement 


y.=-(r--r C ) 


(9) 


This leads to the following relation 


Yi = (Xj - .V, )x + (y-j ~ y, )y + ( - c, )z + 


k f , 2 f \-j 


s 2 \ 1 = 2.3 n 


( 10 ) 


If we consider only n = 3 sensors, (4), (7), and (10) can be combined to give 


~a 2 ' 


V 

tv 

1 

Jr* 

tv 

1 

Vi 

(Z 2 ) 

) 

0 

“j 


(-*}-*,) (y 3 -y,) 
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0 
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where <5 is a weighting factor chosen as 10 3 m, t//j = — 1 /2 is a dimensionless 

parameter, d t = ct, d s = £ 2 .v 2 , and ^ = [F, ( x 2 — .v,) 1 ^] 1 is a scaling factor to be described below. 
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Defining the column vector on the left-hand side of (1 1) as g, (he matrix by K, and the remaining 
vector by f, we may rewrite ( 1 1 ) as 


g = Kf 


( 12 ) 


All elements of K and fare in units of meters, and all elements of g are in squared meters. 
This was accomplished by retaining a factor of c in front of the time difference measurements in the 
first two rows (fourth column) of K, by multiplying (7) by the weighting factor <5, and by scaling 
the field, F, , and source strength, . 9 , each by the factor | , i.e., by making the substitutions: F, -» 
£F, , and s -» £, s in (8). (Note that the sites must be numbered in such a way that .t ; > .v, so g is 
not complex; this can always be accomplished since the numbering of sites is arbitrary and because 
the translation and rotation of the xy coordinate system used in the LP method is arbitrary.) 

In general, the K matrix has 3/1—2 rows and 5 columns, where n = 1, 2, 3, .... If there is 
only n — 1 sensors in the network, K degenerates into a row vector and (12) is underdetermined. It 
there are n = 2 sensors, K will have 4 rows and ( 1 2) will still be underdetermined. For n > 3 sites, 
K will have > 7 rows and (12) will be overdetermined. For overdetermined systems, f can be 
retrieved using the least-squares inversion provided in Twomey (1977) 


f = fKK)-'Kg , 03) 

where the tilde represents matrix transposition. The source time of occurrence and source strength 
are determined as: t — djc , s = d / £ , respectively. 

From the foregoing generalities, we now note that ALDF sensors arc intended to trigger 
only on the ground wave of a cloud-to-ground discharge. Hence, the source can be regarded as 
being located at z = 0- In this case, we can remove the third component of f, that is, we consider 
the column vector f = col(x, y, d t , d % ), and we remove the third column of K. We then regard the 
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expression in (12) as a (3n-2) by 4 system of linear equations. In this case, n = 2 sensors 
generates a (4x4) K matrix so that (12) is a determined system with direct solution f = K~‘g. 
Hence, the LP method can be used by an experimenter that has only two sensors, each measuring 
bearing, arrival time, and field amplitude. In this case, source location (.v,y), time of occurrence, t, 
and source strength, 5 , can be retrieved. If the two sensors do not provide field amplitude 
information, the experimenter can still retrieve the flash location and time of occurrence, i.e., (12) 
becomes a (3x3) system of linear equations, and f = col(x,y,d t ). 

If § is unity, the row vectors of K involving sin <p, and cos 0, appear numerically small i.e.. 
like a zero vector, relative to the other row vectors of K, and the matrix is ill-conditioned for many 
source locations when only three ALDF sensors are available. 

To avoid unstable inversions associated with an ill-conditioned K-matrix, we have made 
the assignment <5 = 10 3 m. This increases the magnitude of the small trigonometric components of 
K and effectively filters small eigenvalues; see section 3.4 below and Appendix A for additional 
details regarding the value of 8. Other, more sophisticated means of filtering small eigenvalues by 
adding external physical constraints to the solution process are discussed in Twomey (1977). 
Chapter 6. 

3. Simulated Tests of the LP Method 
3.1 Overview 

In this and all tests to follow, we do not consider the field amplitude data, F t . We also 
assume that all sources and sensors arc located on the surface of a spherical Earth. By selecting a 
known source latitude/longitude location, we generate the true arrival times and bearings to each 
sensor. Simulated measurements are generated by adding errors to the computed arrival times and 
bearings. The errors are chosen from a uniform random distribution. 

Next, the simulated measurements are analyzed with the LP method. Since the LP method 
is a planar model, we must establish a convention for mapping source and sensor locations 
(expressed in degrees of latitude and longitude) to locations in the .yy-plane of a standard Cartesian 



coordinate system. We then apply the LP method to solve the problem in the Cartesian system. 
Next, an inverse mapping is used to convert the (x, y) solution back into latitude and longitude 
coordinates on the surface of the Earth. At this point, the latitude/longitude solution can be 
compared with the known source to assess true location error. 

If one assumes a flat Earth and performs the entire simulation within a Cartesian coordinate 
system, the resulting retrieval errors are smaller. This is because one avoids errors due to Earth 
curvature and the numerical errors associated with spherical/Cartesian system mappings. Because 
in any real field experiment the source retrievals are ultimately referenced to the spherical Earth, we 
include the net effects of Earth sphericity in this and all other simulations provided below. 

3.2 Spherical Arrival Time and Bearing 

Figure 2 indicates how to compute the arrival time, and bearing, 0,, for the /' ,h sensor on 
a sphere. The unit vectors pointing from the origin O: to the i lh sensor (M), to the lightning source 
( L ) , and to the North Pole (AO are. respectively, 

r = cos cos A u + cos (p t sin A ( v + sin (p i w 

r = cos <p cos Au + cos <p sin Av + sin (pw (14) 


Using the law of cosines from spherical trigonometry gives the spherical angle. A, : 


A = 


180 

n 


-cos' 


( cos a - cosb j cos c . 
sin A sin c. 


(15) 


where a = cos -l (r^ • r), b j - cos"' (r, v • r. ), c. = cos 1 (r • r) - -^cr , R - radius of Earth, and 
by convention the lightning source activates at / s 0. The angle. A,, varies between 0° (North) and 
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180° (South). Wc correct A, to construct the bearing function. <p> , that varies in the manner: 
0° (East), 90° (North), 180° (West), 270° (South), that is, 

<p. = 90 - A. (North East sources) 

<j p =90 + A i (North West & South West sources) ( 1 6) 

<f>. = 450 - A i (South East sources) . 


3.3 Mappings 

In general, different mappings produce different retrieval errors. We consider two possible 
approaches: Mapping #1 (chosen for its mathematical simplicity), and Mapping #2 (chosen for its 
orthogonality). In Mapping #1. we have 

.v = (A - A, )Rcos(p { 


where (A, <p) is an arbitrary longitude and latitude, respectively. The origin of the Cartesian 
coordinate system has been arbitrarily selected as site / = 1, i.e., the ordered pair (A! , (fi\) is the 
location of site 1 and x( A = A, , (p = (p } ) = 0, y (A = A, , (p = <p, ) = 0. Note that y is measured 

along a great circle, i.e., a longitude belt, but that.? is measured along a latitude belt (which is only 
a great circle if <p t = 0 (the Equator)). 

In the second approach, or Mapping #2, we insist that both x and y are measured along 
great circles. To do this, we consider an orthogonal system (u, v, w) where u is a unit vector 
directed from the center of the Earth to the intersection of the Prime Meridian and Equator, w is a 
unit vector directed from the center of the Earth to the North Pole, and v completes the ordered 
triple in accordance with the right-hand rule, i.e., v = w x u. We then rotate this coordinate 
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system through two Euler angles (A, , <p,) and define the new resultant (“starred”) system as 
(u*, v*, w*). In the starred system, u * is directed from the center of the Earth to site 1 . Mapping 

#2 is then 


x = RX*(X,(p) 
y = R(p*(X,(p) , 


(18) 


where 


A* (A ,(p) = tan 1 


cos (p sin A cos A. — cos (p cos X sin X ] 


cos <p cos A cos A, cos <p ] + cos (p sin A sin A , cos (p , +sin<psin<p [ j 
(p“ (X,(p) = sin ~ ' [sin(pcos(p i — cos (p cos Xcos X x sin(p ] — cos (pstn X sin X ] siiup , ) 


(19) 


Again, one can verify from (18) and (19) that ,r(A - A, , <p - (p ] ) - 0, y(A - A, , (p — <j0 ( ) — 0. 
The arctangent expression in (19) must be appropriately corrected depending on what quadrant 
(North East, North West, South East, South West) the point (A, (p ) is relative to site 1. 

3.4 Simulation 

We first consider three ALDF sites in the Darwin, Australia region that were used as part of 
the MCTEX described in Keenan et al. (1994; 1996). Computer-generated lightning sources were 
spaced 0.02° (~ 2 km) apart across the analysis region. Figure 3(a) shows the spatial distribution 
of the retrieved horizontal source location error, contoured in units of kilometers, when no 
experimental errors are considered and when Mapping #1 is used. The retrieved location errors arc 
within 1 km for regions inside the ALDF network. 

Since no experimental errors have been added to the simulated values of the arrival times 
and bearings, the retrieval errors shown in Figure 3(a) are due solely to Earth curvature and 
numerical truncation error. We originally performed these simulations, and all simulations to 
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follow, assuming a flat Earth. When this was done, all of our methods gave retrieval errors well 
below 2.5 m across the entire analysis region when no measurement errors were involved. This is 
to be expected since our methods arc exact solutions for the plane (the 2.5-m error maximum 
occurred only in the LP method over a limited portion of the analysis region and was an artifact of 
what accuracy level we required of our iterative matrix inversion routine). Hence, compared to the 
negligible errors obtained from the flat Earth simulations, the errors shown in Figure 3(a) are 
effectively due to Earth curvature alone. However, the amount of retrieval error due to Earth 
curvature depends on what cartesian-to-spherical coordinate system mapping is used (e.g.. 
Mapping #1, or Mapping #2). Because one will always be interested in how much retrieval error 
the planar models acquire due to Earth curvature, all simulations below show retrievals first 
without added measurement errors, as in Figure 3(a). 

When experimental errors are included in the simulation, we obtain the retrieved location 
errors given in Figure 3(b). The retrieved errors are mean values obtained from performing 100 
individual retrievals at each trial location. For each of the 100 trials, an arrival time error selected 
from a uniform random distribution (ranging from —300 ns to 300 ns) is added to the arrival time 
value, and a bearing error (ranging from -2° to 2°) is added to the bearing value. In addition, we 
have simulated sensor location errors by purposely entering into the LP method false site locations 
(with an error as great as '/: m); the sensor location errors have remained fixed for all source 
analyses. As expected, the addition of experimental errors increases location retrieval errors, but 
the retrieved errors arc still within 10 km for a large portion of the analysis region. Roughly 
speaking, i.e., not accounting for Earth curvature errors, truncation errors, or other errors due to 
matrix inversion, a 300-ns timing error multiplied by the speed of light gives only a 90 m error, 
and a 2° error at a range of 300 km is about 10 km. 

When Mapping #2 is used instead of Mapping #1, we obtain the results shown in Figure 
4(a) and (b). As before, no experimental errors have been added to the sensor positions, arrival 
times, and bearings in the results of Figure 4(a), but the results in Figure 4(b) include these errors. 
The results in Figure 4(a) appear somewhat better than those in Figure 3(a), but the results in 
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Figures 3(b) and 4(b) are similar since the simulated experimental errors tend to mask differences 
between Mapping #1 and Mapping #2. 

Figures 5(a)-(b) and 6(a)— (b) show all of the same type of analyses just described, but for 
the case of four ALDF sensors. The additional sensor clearly helps reduce retrieval error. In 
addition, Mapping #2 produces smaller retrieval errors than Mapping #1. 

When four sensors are available, bearing data are not needed to obtain lightning location 
retrievals (see (11)). By removing the bearing data from the four-station MCTEX region 
simulations, that is assigning <5 = 0, and by applying Mapping #2 we obtained virtually the same 
results as those presented in Figure 6(a) and (b). This is because the weighting factor 5 introduced 
into (1 1) to generate Figure 6(a) and 6(b) is relatively small (i.e., 10’) so that bearing data has little 
influence on the final solution. We also find little change in the solutions for intermediate values, <5 
= 10, I0 : . However, as we increase 5 from 10' to 10\ 10\ and 10 b the retrieval errors increase 
(see Appendix A for more details). 

For a three-sensor ALDF network, bearing data plays a more profound role in helping to 
constrain the solution space. It 5 = 0 (no bearing data used), there would be fewer constraint 
equations than unknowns, and one would not be able to obtain a solution using the LP formalism. 
(Note: a different formalism to be described in section 4 below can be used to find solutions over a 
substantial portion of the analysis region using just three arrival time sensors.) As noted above, 
when S = 1 the K matrix is ill-conditioned for many source locations, and the computer time 
required to invert K is excessive. When 5 = 10, 10 2 , 10\ or 10 4 there is no problem inverting K 
and there is no appreciable change in retrieval error. When 5 = 10 5 the retrieval errors begin to 
increase slightly. Three-station LP simulations for different values of 5 are also provided in 
Appendix A. 

For comparison, we also provide error results (Fig. 7(a) and (b)) for the three sites used in 
the Tropical Ocean Global Atmosphere Coupled Ocean-Atmosphere Response Experiment (TOGA 
COARE) described in Orville et al. (1997) and Peterson et al. (1986). This experiment employed a 
larger sensor baseline than that used in MCI EX, and our simulated tests cover an analysis region 
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1 8°x 1 8° in latitude and longitude. The sources in this simulation where placed 0.05° apart and 
Mapping #2 was used; as with the three-sensor MCTEX study, the retrieval errors for Mapping #2 
differ little from those errors obtained using Mapping #1. Overall, the 2° bearing error and the 
effects of Earth curvature make it difficult to obtain errors below 10 km for distant sources. 

4. Quadratic Planar (QP) Method 

In this section we assume that only three arrival time measurements are available from the 
ALDF network. Hence, the methods of section 2 cannot be applied, but some insight about the 
source location can still be obtained. Once again, considering only the radio emission from the 
lowest part of the return stroke, we take z = 0. We assume that sensor / = 1 is at the origin ot a 
rectangular Cartesian coordinate system, i . e . . .v ( = v ( = c, =0, and we specify the convention 
/, =0. This leads to one nonlinear equation and two linear equations all in the three unknowns 

(•*. v, r) 


x- + V* = r- 


q'i = j(rf- - c 2 tf ) = XjX + y,y + ct { r ; i = 2, 3 


( 20 ) 


The equations in (20) were derived from the transit equation in (1) by means similar to that 
discussed prior to (2) of section 2. We have removed the source activation time, t, with the 
relation t = -r/c. Since r, a 0, it is consistent that t < 0 given that r > 0. Each sensor is a distance 

r. = ( x f + y ( 2 ) ,/2 from the origin. 

Fundamentally, these equations define three circles each with radius c ( /• - 1) and centers at 
(jc where i = 1, 2, 3. We will see below that the source is located where the three circles 
intersect. We will also find that certain source locations produce arrival time data that can be 
described geometrically by two possible sets of three circles. Each set of three circles define a 
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unique intersection point in the .vv-plane, thereby leading to a fundamental ambiguity in source 
location retrieval. 

Geometric intersections of the circular curves described above are obtained by solving the 
system of equations in (20). To solve the system, we first subtract the terms c/,r from each side of 
the linear equation set 


4,-1 = 4/ “ ct i r = x jX + VjV ; i = 2.3 . 
Identifying the vectors, q = col(r/, , c/,), r = col(.v, y), we may write 


( 21 ) 


q =Qr 


( 22 ) 


where the Q-matrix and it’s inverse are given by 


X 2 

r ‘ 

ri 

J**-. 

D 

1 

1 

' >3 

-y 2 " 

_ x i 

■V 

• v — 

-v ; y 3 -y 2 .v 3 

_-*3 

* V 2 


(23) 


From (22) and (23) and our discussion preceding (2 1 ) we have the relations 


x( r) = (\y/, - y 2 q 2 )/(x 1 y J - y 2 .r 3 ) 

}'(r) = (x 2 q 2 - -ty/, )/(x 2 y i - y 2 x 3 ) (24) 


t(r) = -r/c . 
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The x and y variables are written as functions of r in (24) since the components of q depend on r as 
given in (21). Substituting the first two equations of (24) into the first (nonlinear) equation of (20) 
and carrying out the algebra leads to an equation quadratic in r alone 


Ar 2 + Br + C = 0 , 


(25) 


where 


A = c 2 [r 2 / 2 - 2(.v,.t 3 + y,y, )f 2 / 3 + r, 2 / 2 ] - (.t, y 3 - v 2 .v 3 )' 

B = 2 <^ <■}■>! 2 y 2 -^* 3 ^ 3 ^ 2 ) — ^2 ^ 3^3 ] 

C = r 3 2 <7 2 2 - 2(.v 2 a 3 + y,y 3 )q\q\ + r}q’ y 2 . 


(26) 


Hence, the lightning source range, r, is the nonnegative real root obtained from the formal (two 


root) solution 


-g + y/i? 2 -4AC 

o a 

( 21 ) 

-fi-3/g 2 -4AC 

2A 


Values of r = 0 correspond to a direct lightning strike of sensor i = 1, which we ignore. Note from 
(26) that the numerical value of the coefficients (A, B, Q are obtained from the sensor locations 
and excitation times, that is, on the six variables: {x 1 ,y l ,t 1 , x v y v t^-, the variables, q', are 
obtained from the expressions y(r 2 - c 2 / 2 ) as given in the last two equations of (20). After these 
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data arc used to compute r, (24) is used to find the lightning location (.v(r), y(r)), and time of 
occurrence, t(r). 

5. Simulated Tests of the QP Method 

By placing computer-generated lightning sources 0.02° apart in latitude and longitude 
across the analysis area (see section 3), we have determined the horizontal location error resulting 
from each root in (27). To facilitate comparisons with simulated tests ot the LP method, sensor 
position and arrival time errors used here are as described in section 3.4 and 100 trials at each 
source location are once again used to generate mean retrieval location errors. Mapping #1 given in 
(17) is employed. 

Figure 8 clarifies what root provides a smaller retrieval error. The shaded regions are where 
r + provides a better retrieval than r (the unshaded region is where r provides a better retrieval than 
rj. Interestingly, the dividing lines of these regions are defined by the sensor baselines. 

When we pick only the root that provides the best retrieval and plot the associated error 
result over the analysis region, we obtain the result given in Figure 9(a). When a 300-ns uniform 
random error is added to the computer-generated arrival times, we obtain the mean horizontal 
distance errors given in Figure 9(b). Considering that only three sensors are involved, retrieval 
errors are quite good; a large region ot errors below 1 km is evident. Distant sources, or sources 
located near the outer sensor baselines are more difficult to accurately retrieve. By outer sensor 
baseline we mean: any position along a great circle passing through two sensors except those 
points on the great circle located between the two sensors. A comparison between the three-station 
LP results in Figure 3(b) shows that the QP method provides better results over most of the 
analysis region. This is due in part to the large 2° bearing errors implemented in the LP simulation 
and the fact that the three-station LP method depends on bearing data to obtain a solution. The 
four-station LP method (Fig. 5(b)) outperforms the QP method. 

Since bearing data can aid in determining which root, r + or r_ , is correct (see section 6. 1 
below on solution ambiguity) and since the QP method gives generally better results than the three- 
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sensor LP method, it is evidently better to use the QP method than the three-station LP method 
even when bearing data are available. This conclusion is based, ot course, on an assumed bearing 
error of 2". 

6. Examination of QP Method Roots 

When applying the QP method to actual arrival time data, one picks the solution associated 
with a nonnegative real root, i.e., the source range, r, must be nonnegative and real. A detailed 
discussion of root results is provided below. 

6.1 Unequal Nonnegative Real Roots (Ambiguities) 

In the discussion of three sensor networks by Hollc and Lopez ( 1993), pp. 8 and 11. the 
lightning source location (a. y) is described in terms of the intersection of two hyperbola branches; 
each branch is defined by two sensors. For some lightning source locations, the hyperbola 
branches intersect at two locations (see for instance. Fig. 6, p. 11 of Holle and Lopez, 1993). This 
amounts to a fundamental ambiguity in location retrieval and the authors correctly assert that the 
ambiguity can be removed by adding a fourth (properly positioned) sensor. Ambiguities are 
described in our formalism by the intersection of circles as indicated above following (20). To fully 
appreciate the hyperbolic and circular geometrical viewpoints, it is important to recognize that the 
two intersection points defined by two sets of three circles are identical to the double intersections 
obtained from the two hyperbola branches mentioned above. In other words, these two widely 
different geometrical viewpoints produce identical results, as they must. 

Without reference to the geometry of hyperbolic or circular intersections, our algebraic 
formalism immediately defines all ambiguous cases. An ambiguity will exist whenever two 
unequal nonnegative real roots result from (27). In order to determine what lightning source 
locations produce these “ambiguity regions”, we have kept a record of the root results in the 
numerical experiments described in section 5 above. For the case of no simulated experimental 
errors, the source locations that resulted in two unequal nonnegative real roots are indicated by the 



shaded regions in Figure 10; see section 6.2 below for minor corrections to the ambiguity regions. 
In general, a different network geometry would produce different results. 

Strictly speaking, since two distinct sources can produce identical arrival time difference 
information, there is no means of discriminating which source location is correct unless some 
additional information is supplied to the retrieval process. In effect, the solution is fundamentally 
nonunique. (Similar comments about the retrieval ot charge from ground-based held measurements 
have been made in Koshak and Krider ( 1994). In that problem, a point charge Q„ , and a sphere ot 
radius a with total charge Q, produce identical electrostatic fields outside the radius a.) Hence, 
additional measurements (e.g., arrival time, bearing, signal amplitude, radar, acoustical, 
interferometric) must be used to pick the correct root. Bearing data would be the most common 
data to use in root discrimination since it is part of the ALDF data stream. 

Nonetheless, an experimenter might be tempted to compare the shaded regions in Figures 8 
and 10 in order to determine which ot the two unequal nonnegative real roots produce the true 
source location. However, one must remember that Figure 8 does not pose a real physical 
constraint to an unknown source since it only provides information if the source location is 
already knowtv, obviously in a real field experiment the source location is not yet known. 

Additional rigor clarifies the immutability of the ambiguous case. Note that there are three 
ambiguity regions in Figure 10, and call the regions |/?|, R i, /?}}. Similarly, there are three regions 
{ , P 2 , P 3 } in Figure 8 that are wholly contained within each of the respective ambiguity regions. 
If one subtracts the respective regions { P, , P, , P y } from the respective regions {R f ,R 2 ,R } }, one 
obtains the three regions {N^N^N^s {/?, - P { , R 2 - P, , /? 3 - P 3 }. When two unequal 

nonnegative real roots are obtained, we find that r + occurs in P , and r_ occurs in ;V,, that is, the pair 
of roots produce solutions in (P r Nf), (P 2 ,N 2 ), or (P^,N^). From Figure 8, region P , is the 

region where r ¥ is correct, N-, is the region where r_ is correct, and both P , and N, are subregions of 
R h In other words, each solution is possibly correct so that comparisons between Figure 8 and 
Figure 10 serve no help in determining the correct root. 



Nonetheless, r t is correct for most of the ambiguity region shown in Figure 10. Hence, an 
experimenter who obtains two nonnegative real roots, but who does not have ancillary data sets 
such as radar, magnetic bearing, etc., has a better chances of getting the correct solution if he/she 
selects the root r, . 

6.2 Equal Nonnegative Real Roots 

In this section we are interested in identifying what source locations produce two equal 
nonnegative real roots. Note that this condition is satisfied when the discriminant, 
B^-AAC, in (27) is zero, that is, the two equivalent roots correspond to a unique (unambiguous) 
solution [,v(r), y(r), t(r)] where r + = r_ = r. In Appendix B, we show for arbitrary noncollinear 
network geometries that the discriminant function is zero only along the outer sensor baselines. 
Therefore, the ambiguity regions shown in Figure 10 are technically not ambiguous along these 
linear domains. 

6.3 Complex Roots 

Complex roots occur whenever the discriminant in (27) becomes negative. Figure 1 1 
shows how the discriminant varies for different source locations across the analysis region. As we 
have already shown in section 6.2 and Appendix B, the discriminant is zero for sources located 
along the outer sensor baselines. Figure 1 1 shows additionally that the discriminant is a relative 
minimum at the outer sensor baselines. 

From the simulation in section 5 (with 300-ns arrival time errors, and 100 trials per test 
location) we have tallied the fraction of trials at each location that produce complex roots. Figure 12 
shows that there are no complex roots over most of the analysis region except when the sources are 
near the outer sensor baselines. These regions (or “spokes”) appear to diverge with range from the 
sensors and as many as 40-60% of the sources are complex within the spokes. Clearly, for 
sources located sufficiently close to the outer sensor baselines, measurement errors are 
occasionally large enough to drive the discriminant negative. Whenever the discriminant is 
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negative, both roots in (27) are complex and no physical solution is obtained. Conversely, 
whenever complex roots are obtained from a set of actual measurements the source is likely to be 
located in one of the spoked regions. 

6.4 Overview of Root Results 

From our discussion so far, we can conclude that any retrieval will produce one of the 
following cases: (a) r+ > 0, r_ > 0, r + r_ ; (b) r + > 0, r_ > 0, r+ = r_ ; (c) /%. <0, r_ > 0: or (d) r+ 
and r_ complex. Case (a) corresponds to a source that is located inside one of the ambiguity 
regions, case (b) corresponds to a source located on an outer sensor baseline, case (c) corresponds 
to a source that is not located in any of the ambiguity regions or along any outer sensor baseline, 
and case (d) corresponds to a source located on or near any outer sensor baseline when 
measurement errors are sufficient to drive the discriminant negative. 

Note that we do not include the case r + > 0, r_ < 0 since if >\ > 0. the source must lie in one 
of the ambiguity regions implying that v _ would be nonnegative (i.e., a contradiction). We also 
disregard the case that both roots are negative since a physical source must lie a nonnegative 
distance from sensor i = 1 . 

7. Sample Storm Analyses 

We have applied the LP and QP algorithms (with Mapping #2) to retrieve several thousand 
cloud-to-ground flashes that occurred over the MCTEX analysis region during November 28 (Day 
332) and November 29 (Day 333), 1995. To demonstrate the internal consistency between various 
forms of the algorithms and to clarify to what degree the algorithms differ, we have analyzed only 
those flashes that produced a valid excitation at all four of the ALDF sites. In this way, each flash 
can be analyzed by each selected algorithm. 

The LP two-station algorithm, or “LP2” algorithm ingests the arrival time and bearing data 
from sites 1 and 2, and ignores the data from sites 3 and 4. The LP3 algorithm ingests the arrival 
time and bearing data from sites 1 , 2, and 3, and ignores site 4. The QP algorithm ingests only 
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arrival time data from sites 1, 2, and 3, and ignores site 4. Finally, the LP4 algorithm ingests 
arrival time and bearing data from all lour sensors. Hence, the results derived from LP4 arc 
considered to be the best retrieval results against which the three remaining algorithms are 
compared. For example, the discrepancy between LP2 and LP4 results indicates how appropriate it 
might be (in future storm analyses) to use LP2 when only two sites are available. 

Figure 13 shows the ground flashes derived from Day 332 for each algorithm, and Figure 
14 shows the results derived from Day 333. Note from the LP2 results provided in Figures 13(a) 
and 14(a) that there are obviously location errors along the baseline between site 1 and 2, in 
agreement with the general error results provided in section 3. The cluster ot flashes to the south 
on Day 332 shown in Figure 13(d) are reasonably well retrieved by LP2 (Figure 13(a)), but are 
somewhat smeared out and biased closer to the network. Similarly, the myriad ot flash clusters to 
the south on Day 333 shown in Figure 14(d) are reasonably well retrieved by LP2 (Figure 14(a)), 
but again there is positional smearing. 

The results of LP3 (Figures 13(b) and 14(b)) dramatically reduce the baseline errors and 
smearing that was associated with the LP2 results. However, there are still some positional 
adjustments between the LP3 and LP4 results. 

Finally, the results of the QP algorithm are shown in Figures 13(c) and 14(c). In all of the 
QP results shown here, root ambiguities were resolved using the LP3 source solutions. Generally, 
these results agree favorably with the LP4 results and are perhaps even better correlated to the LP4 
results than arc the LP3 results, for some source locations. Of course, some hyperbolic-shaped 
artifacts are evident in the QP results in agreement with the error results provided in section 5, and 
several complex root solutions had to be removed. Because the LP4 and QP methods are 
significantly different mathematical approaches, it is particularly encouraging to see such good 


agreement between them. 
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8. Summary 

In this writing we have derived, tested, and applied two basic methods for retrieving the 
location and time of occurrence of lightning ground strikes from a network of ALDF sensors. We 
developed the methods so that they could be used in future validation studies of NASA space- 
borne lightning sensors: the OTD and the LIS. The validation requirements of OTD and LIS have 
allowed us to avoid the expenditures associated with using the proprietary nonlinear X' 
minimization algorithm mentioned in Cummins et al. (1993, 1995, 1998). The development of our 
algorithms has also expanded upon and clarified the theoretical aspects ot ALDF data inversions. 

The first approach introduced in this writing, or LP method, assumes that arrival time, 
bearing, and field amplitude measurements are all available from the network. As provided in (12), 
these measurements are collected into one coherent linear system of equations that is solved by 
straightforward inversion. Because ot the general torm ot (12), we have clarified what solution 
options one has if only a subset of the measurements are available and/or if one or more sensors 
could not trigger on an event. In the extreme case of having only three arrival time measurements 
from the ALDF network (as might be the case due to sensor subsystem failure), we have 
introduced the QP method. 

The planar methods express the source locations directly in terms of the measurements. The 
solutions are concise, require little computer time, and afford the user with specific physical 
insights about the retrieval problem (e.g., relative importance/effects of timing/bearing data on final 
solution, regions of ambiguity, source locations producing complex roots). This starkly contrasts 
the approach of nonlinear x 2 minimization wherein solutions are found by a computer scorch of 
the optimum lightning source parameters (latitude, longitude, and time of occurrence). 

Moreover, the LP and QP methods introduced here offer the authors and other researchers 
a means to intensively analyze and compare, first hand, lightning radio source locations with 
OTD/LIS low-Eaith orbit lightning detections. In the future, we intend to apply these methods to 
analyze a wide range of thunderstorms, to continue intercomparing the methods, and to relate the 
results to OTD and LIS and other independent data sets such as: radar. Lightning Detection and 
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Ranging (LDAR), and the NLDN. The first author will also improve some of the matrix methods 
presented here to directly account for Earth sphericity; more elegant oblate spheroidal models are 
also under consideration. 
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Appendix A — Weighting of Bearing Data 

It is of interest to deteimine to what extent bearing data is actually being used to constrain 
the lightning source location in the LP method as a function of the weighting factor <5 introduced 

into the linear system of equations provided in (1 1). Insight is gained by considering the case of 
three sensors and ignoring field measurements, F . Then (11) reduces to 


’ a 2 ' 


’(.v ; - .r,) (y 2 - y, ) c(7, -t z ) 



(.y,-.v,) (y, - y, ) c(r l -/,) 

<5A 

= 

5 .sin 0 | -8 cos <p\ 0 

<5/T 


8 sill 02 ~8 COS (p 2 0 

M. 

! 

0 sin 03 - 8 cos 03 0 


This system, which can be written in the standard notation: g = Kf, has the least squares solution 
f = (KK)“'Kg. We are interested in the explicit functional dependence of x and y on the arrival 
time. /, , and bearing, 0, , data, and on the weighting factor, S. Because this is a very involved 
hand calculation, we utilize a computer-aided symbolic manipulator to arrive at the following torm 


x = 


li^d 1 + / j , 


/liQt-i + /q « 3 + ^(hjS 2 +hj +3 )Pj_i 

;=•» 


(A-2) 


The ten functions, h J (j = 0, ... , 9) depend, in general, on the arrival time and bearing data. The 
variables (a, , «,) depend only on arrival time data and network geometry, and the variables (J3, , 

/J 2 , (i } ) depend only on bearing data and network geometry as given in (3) and ( 6 ), respectively. 
A similar form holds for y. The coefficient in front of the square brackets in (A-2) docs not 
preferentially weight the a ' s or /Ts so it is of no concern in this discussion. However, the 
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coefficients (h j +S z h^ 3 ) do weight the p's but not the a's (i.e., these coefficients weight the 
bearing data, but not the arrival time data). This leads to the final results 


5 = 0 => no solution 

lim x = y^-(/i 4 /J| +h 5 p2 / = !, 2, 3. 

S—>°o '0 


(A-3) 


The first result is true since when <5=0 there are two equations in three unknowns in (A-l). The 
second result is true because the only arrival time dependence associated with each function 
( Ijq , hj , hj , h t ) is a factor (/; + /*) . When the ratios ^//zyare taken this factor cancels out and we 

are left with a function. A, that depends only on bearing data. 

In summary, for a positive finite value of 5, both arrival time and bearing data are utilized. 
However, as 5 is increased from zero, bearing data eventually becomes more heavily weighted 
over arrival time data until, for a sufficiently large value of 5 , only bearing data is being used to 
determine the source location. For 5=1, we have difficulty inverting K for many source locations. 
We have also performed retrievals, in the presence of measurement errors, for the values: 5 = 10, 

10 2 , 10\ 10\ 10 5 , and 10 6 . Figure A. I shows the result for <5 = 10, 10\ and 10 5 . There is not 

much change in the solution from 10 to 10\ but the dominance of bearing data constraints at 10 5 
begins to reduce the quality of the solution (i.e., a 2° random bearing error can create a substantial 
location error if the source range is sufficiently large). 

We have performed the same type of computer-aided symbolic manipulation to determine 
explicite forms when a four-sensor network is used (i.e., one more arrival time equation and one 
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more bearing equation is added to the system in A- 1 so that K becomes a 7x3 matrix). In this case, 
the form of x (as well as y) is 






(A-4) 


Most of the 17 functions, k , j = 0 16 depend on both arrival time and bearing data. 

However, (/r-,,/r 6 ,A' 7 , A' s ) depend on arrival time data, but do not depend on bearing data. Wv 
obtain the final limiting conditions 

lim .v = 7 — (kfrCCi + k-jCt j + k^ttj ) = r(tj ) 

5 -° ' 2 “ (A-5) 

lint x = T—[kqP) + k\QP^ + £| + k^Pi) = ) ; i = \,2,3 . 


Hence, we swing from a solution governed only by timing data (5 = 0) to one governed only by 
bearing data (6 = °°). The solutions for 5= 10, 10\ and 10 s are shown in Figure A. 2. The 
bearing data significantly worsens the solution when weighted heavily ( 8 = 10 5 ). 
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Appendix B — Locations Where QP Method Discriminant Function Vanishes 

We investigate more rigorously the zeroes of the discriminant function A = (B 2 -4AQ of 
(27) in the QP method. The computer plots of this function gave some interesting results near the 
outer baselines of the sensors, i.e., there appears to be minima there (see Fig. 11). 

In the following formal approach, we algebraically reduce the discriminant into the product 
of three factors. Each factor is then shown to vanish along a specific outer sensor baseline. Our 
results apply to arbitrary network geometries. Because a zero discriminant implies that two 
nonnegative, real, and equal roots are obtained, a unique (unambiguous) solution (-t(r). y(r), f(r)) 
is obtained on the outer sensor baselines where r = -B/(2A) = r + = r_. 

Using the forms in (26) for .A, B , and C, the discriminant can be written as 


A = 



r 2 r 3 )w ~~ ( l 2 e h t 2 t ?> +( l3~ , 2 ) + > 


( B- 1 ) 


where 

p = iy r, = a,a 3 + y 2 y 3 

£ = -del Q = -(as v 3 - vsjc 3 ) (B-2) 

c c " • 

0 = r 3 -q 2 - ~2pq 2 q i + r±q 3 ~ . 


To simplify some of the algebra without losing generality, we rotate the x and y axes so that 
y 2 = 0. Further reduction of (B- 1 ) leads to 


A = [ x \yl{ x l -c 2 /|)][r 3 2 -c 2 /|][(x| - c 2 t\) + 2(c 2 / 2 / 3 - * 2 .v 3 ) + (r 2 . (B-3) 
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Figure B.l considers the second factor, (r, 2 -c 2 / 2 ), in (B-3). For a source located on the solid line 
with r > r, site three is excited at 


h 



(B-4) 


The factor becomes 



-c-'i) 







(B-5) 


Proceeding in a similar fashion, the factor for a source on the dashed line with 0 < r < r, is the 
nonzero result, 4r(r 3 -r), and the factor for a source on the dotted line (with 0< r < r, ) or a 
source on the thick line (with r>r,) is zero. Hence, the factor is zero along the line running 
through the sensors (including the sensor locations themselves, but excluding the line segment 
between the sensors). This is what we refer to as the “outer sensor baselines”. Similar comments 
can be made regarding sites 1 and 2 when the factor (.v; - c 2 f; ) is considered. 

Evidently the third factor in (B-3) corresponds to the line running through sites 2 and 3. To 
prove this, we consider the geometry provided in Figure B.2. For a source on the solid line with d 
> 0 we have 


D + d-r 


d-r 

h = 


(B-6) 


c 



Substituting these expressions into the third factor, and noting that D 2 - (x 2 - .t, ) 2 + yj - 
.v 2 - 2.r,.v, + r 2 , we obtain 

xj -(D + d - r) 2 + r 2 - (d - r)~ + 2( D + d - r)(d - r) - 2.V-..V, 

= (x; ~2 x 2 x^ +r 2 )-[(D + d-r)-(d-r)Y (B-7) 

= (x; -2.v : .r 3 +r~)-D 2 =0 . 


For a source located on the dashed line (but not at site 2 or site 3) the factor reduces to the nonzero 


result: 4a) 2 (O y (.t, -.r,)' +.v 2 , where the constant factors (o) 2 ,co y ) obey the constraints: 
G), + ft), = 1 . <0-, > 0, ft), > 0 . Finally, for a source on the dotted line and a distance / > 0 from site 
2, we have 


t. 


I - r 


c 


( B-Sl 


D + l-r 


These expressions have the same form as those in (B-6) but are interchanged. When substituted 
into the third factor of (B-3), the factor reduces to zero as in (B-7). This completes the proot 
showing that the discriminant function vanishes along the outer sensor baselines. 



33 


Figure 1. Geometry associated with the LP method. 

Figure 2. Spherical trigonometry used for determining arrival time and bearing. 

Figure 3. Lightning location retrieval errors for a three-station network using the LP method: (a) 
no measurement errors, and (b) with the following random measurement errors: 0.5-m sensor 
location error, 300-ns timing error, 2° bearing error. Mapping #1 is used. Contours are in units ot 
kilometers. The analysis region shown is 6° in latitude (667 km) by 6° in longitude (about 651 
km), and this is where the Maritime Continent Thunderstorm Experiment (MCTEX) was 
conducted. 

Figure 4. Same as in Figure 3 except that Mapping #2 is used. 

Figure 5. Same as in Figure 3 except that tour sensors arc used. 

Figure 6. Same as in Figure 3 except that four sensors and Mapping #2 are used. 

Figure 7. Same as in Figure 3 except that this is for the Tropical Ocean Global Atmosphere 
Coupled Ocean-Atmosphere Response Experiment (TOGA-COARE) analysis region, and mapping 
#2 is used. The analysis region is 18" in latitude (2002 km) by 18° in longitude (about 1996 km). 
Shading, rather than contouring, is used to clarify the nonmonotonic distribution of retrieval 
errors. 

Figure 8. Comparison of roots in the QP method. The shaded regions indicate where root r ; 

produces better retrieval results than root r . The unshaded regions indicate where root r_ produces 

better retrieval results than root r + . No measurement errors have been added to the simulated data. 
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Figure 9. Retrieval errors from the optimum root in the QP method: (a) no measurement errors, 
and (b) with the following random measurement errors: 0.5-m sensor location error, 300-ns timing 
error. Mapping #1 is used. Contours are in units of kilometers. 

Figure 10. Shaded regions are the QP method “ambiguity regions’’ that indicate what lightning 
source locations result in two unequal nonnegative real roots. No errors were added to the 
simulated arrival times. 

Figure 11. Plot of the (scaled) discriminant of the QP method. Expressions for A, B. and C in the 
discriminant, B 2 -AAC, are given in (26) of the manuscript. 

Figure 12. Fraction of 100 simulated sources at each location that produce complex roots using 
the QP method. A 0.5-m sensor location error, a 300-ns timing error, and Mapping #1 was used. 

Figure 13. Day 332 ground flash retrievals using the following algorithms (see algorithm 
descriptions in text): (a) LP2. (b) LP3, (c) QP, and (d) LP4. Only those source retrievals falling 
inside the above MCTEX region are shown. Of the four algorithms, the results of LP4 are 
considered the best estimate of the true lightning locations. 

Figure 14. Day 333 ground flash retrievals using the following algorithms (see algorithm 
descriptions in text): (a) LP2, (b) LP3, (c) QP, and (d) LP4. Only those source retrievals falling 
inside the above MCTEX region are shown. Of the four algorithms, the results of LP4 are 
considered the best estimate of the true lightning locations. 


Figure A.l Location retrieval errors for a three-station network using the LP method when (a) 8 = 
10 (b) 8 = 10 3 and (c) 8= 10 s . The same simulated measurement errors discussed in section 3 are 
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used: 0.5-m sensor location error, 300-ns timing error, 2° bearing error. Mapping #2 is used. The 
contours are in units of kilometers. 

Figure A.2 Same as in Figure A. 1, but for a four-station network. 

Figure B.l. Geometry for analyzing the second factor in the discriminant function. 

Figure B.2. Geometry for analyzing the third factor in the discriminant function. 
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